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1 Introduction 

For EEG- and MEG-source reconstruction so-called 
equivalent dipole models are commonly used. All 
equivalent dipole algorithms need a-priori 
knowledge of the number and class of sources 
involved in the brain activity to be reconstructed. 
Anatomical a-priori knowledge of the source 
positions (and orientations) cannot be incorporated 
into these equivalent dipole methods directly. 

During the last years an increasing interest in current 
density reconstruction algorithms can be seen. All 
these methods have in common, that elementary 
dipoles are distributed on given positions, e.g. on 
regular grids inside the head or on the gray matter 
layer. The calculation of the strengths and 
orientations of these dipoles generally leads to a 
highly under-determined system of equations 
(number of unknown dipole components is much 
larger than the number of sensors). It can be solved 
by additional, appropriate constraints only. By 
applying these (e.g. minimum norm) constraints, 
mathematically unique solutions can be achieved. 
The additional, constraining model term has to be 
weighted against the data term, which is done by a 
so-called regularization parameter. 

Current density methods are commonly applied to 
single time-point evaluations, but do not need the 
number of sources as additional input. Also no 
assumptions about the shape or size of an activated 
area are made, whereas single equivalent dipoles are 
based on the implicit assumption, that patches of 
activated gray matter are well represented by their 
centers of gravity and their mean surface normals. 
Another hope is that one can estimate the extents of 
the active brain areas with current density recon¬ 
structions as well. 

Due to the underdetermined character of the problem 
and the blurring properties of the minimum norm 
least squares constraint, this class of source 
reconstruction algorithms yields smooth current 
distributions in contrast to the rather over-focussing 
equivalent dipole methods. 

2 Methods 

For the least squares approaches the dipole 
components are calculated for each time slice by 
solving a linear system of equations using matrix 


inversion or Singular Value Decomposition. If 
another norm is used instead of the L 2 -norm of the 
least squares methods, e.g. the L r norm, a system of 
nonlinear equations with the same number of 
unknowns has to be solved. Therefore, the latter 
methods need much more computational efforts. 
However, the L r norm approaches promise much 
more focal solutions and a more robust behavior 
against outliers in the measured data. 

2.1 Generalized minimum norm formulation 

In a generalized formulation the data term D, that 
accounts for deviations between measured data M 
(column vector of the s sensors) and calculated 
sensor signals S = L*i (with the lead-field matrix L 
(s*(c*d), comprising dipole positions and volume 
conductor properties, and the column vector j of the 
c*d current components), can be written as: 

D=HE*(M : S)llm P _ _ (1) 

The matrix W is an ^-weighting matrix for the s 
sensors, m denotes the norm used for the weighted 
differences and p is the power applied to the data 
term. 

The minimum norm model term T contains the 
weighting matrix C ((c *d) 2 ) for the current 
components, n describes the norm used, and q is the 
power applied to the model term: 

T = HQ*j]ln q ( 2 ) 

Both terms have to be minimized simultaneously 
and are weighted by the positive regularization 
parameter A: 

minimize A= D + A*T = 

/\W*(M - S)//m P + A*llQ*jJln q (3) 

2.2 Minimum norm least squares 

The commonly used minimum norm least squares 
case (m=n=p=q=2) leads to Frobenius norms for 
both data and model terms. The Frobenius norm of a 
it *v matrix A is given by: 
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With non-singular weighting matrices W and C the 
solution of Eq. 3 for the best fit current components 
is in this special case is given by: 

minimize A = //W*(M - L*j)H 2 + ^7/£*i/// (5) 
I = (Q^Q)' 1 *L t *W t *K' *E*M (6) 



with x = E%*(Q T *Q% t *K T + A* 1 (7) 

With no detector and no current component 
weighting (W=l and Q=l) the standard solution for 
the minimum nom least squares case is easily 
obtained: 

minimize A = //M- L*jJl 2 2 + ^*////// (8) 

i = l T *(l% T + h*!)' 1 *M (9) 

This standard solution is known to lead to 
superficial source distributions, since small currents 
close to the detectors can produce fields of similar 
strengths as larger currents at greater depths, which 
will thus be suppressed by the unweighted minimum 
norm term in Eq. 8. To compensate for the undesired 
depth dependency of this approach, the currents can 
be weighted in order to account for the lower gains 
of deeper dipole components (lead-field norma¬ 
lization). 

minimize A = //M - Z,*//// + (10) 

j = (Q T *Q- 1 % T *(l*iQ T *Q- 1 *l T + /W)' 7 *M(11) 
With a location- or component-wise depth weighting 
the large ((c*d) 2 ) matrix C is diagonal and therefore 
(Q^Q)- 1 can be calculated very easily. The other 
matrix inversion (L*(C T *Q~ 1 *L T + A*T)' 7 is applied 
for a much smaller (s 2 ) matrix only and is thus very 
fast, since the number of sensors s is usually in the 
order of 100. 

2.3 Simulations 

We studied different minimum norm algorithms in 
view of their spatial resolution, regularization, and 
lead-field normalization properties. Two classes of 
minimum norm solutions are examined: Linear least 
squares methods and non-linear L r norm approaches 
[1]. Two special cases of linear algorithms, the well 
known Minimum Norm Least Squares and an 
implementation with Laplacian smoothness 
constraints [2], are compared to two nonlinear 
algorithms comprising sparse and standard L r norm 
methods. 

In a Signal-to-Noise-Ratio framework the critical 
choice of the regularization parameter can be auto¬ 
mated with some methods, so that no over- or under¬ 
fitting of the data occurs . Otherwise, this parameter 
has to be adjusted manually, in order to achieve an 
adequate representation of the measured data [3]. 
Compensation methods for the different depth 
dependencies of all approaches by lead-field 
normalization are discussed. Simulations with 
tangentially and radially oriented test dipoles at two 
different noise levels are performed to reveal and 
compare the properties of all approaches. 

Fig.l shows the spherical volume conductor model 
with 81 electrodes and the test-dipole positions used 
for the simulations, whereas Fig. 2 displays the 



Figure 1: Outer shell of the spherical volume 
conductor model (81 electrodes, extended 10/20 
system) and 19 tangential and radial test-dipole 
positions used for the simulations (frontal view). 



Figure 2: Relative Signal-to-Noise-Ratio as function 
of the relative eccentricity of radial and tangential 
test-dipoles normalized to the most superficial 
radial dipole field. 

relative SNRs for tangentially and radially oriented 
dipoles at different depths (eccentricities). 

All current density reconstructions are performed 
with support points on a regular 3D-grid with 7mm 
grid spacing, clipped by the innermost volume 
conductor sphere (84mm radius). This leads to 7497 
points and thus to a lead-field matrix of size 
81*22491. No orientational constraints are used. 

All computations were performed using the Curry 
V4 (Neuroscan) software package. 

3 Results 

Fig. 3 displays the results of three different methods 
at seven selected depths for radially oriented dipolar 
sources. The most superficial source position (99%) 
leads to an SNR of 10, the deepest position (-24%) 
to an SNR of 3 (compare Fig. 2). 

In order to more quantitatively evaluate the locali¬ 
zation and spatial resolution power of the different 
approaches, weighted centers of thresholded distri¬ 
butions (at 50% of the maximum current) are 
calculated and shown in Fig. 4 and 5 for all depths. 
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Figure 3: Different three-dimensional current density 
reconstruction methods for radial test-dipoles at 
different depths with optimal regularization and lead- 
field depth normalization for realistic noise levels. 

The standard MNLS approach results in too 
superficial source distributions, whereas LORETA 
in general overestimates the source depth. The non¬ 
linear methods perform best for sources at large 
eccentricities (> 40%). 

The spatial resolution of LORETA is superior to 
MNLS, especially for deeper source positions. 

Finally, the influence of the regularization parameter 
for a radial dipole at 71% relative eccentricity is 
studied for the two linear methods (MNLS and 
LORETA) in Fig. 6. Overfitting the data (relative 
deviation < 1) results in spurious source distribu¬ 
tions, whereas underfitting the data (relative 
deviation >1) yields a dominating model term 
which in both cases leads to a smooth distribution 
around the center of the volume conductor. 




Figure 4: Quantitative evaluation of different linear 
and nonlinear CDR methods with respect to their 
localization power (weighted center / test-dipole 
location). 



Figure 5: Spatial resolution of the linear CDR 
methods (MNLS, LORETA) for radial and tangential 
test-dipoles. 

4 Discussion 

Linear L 2 -norm reconstruction algorithms yield 
smooth current distributions in contrast to the non¬ 
linear L r norm approaches that result in more focal 
solutions. 

The partial lead-field compensation gives the best 
results for the standard Minimum Norm Least 
Squares and Laplacian smoothness algorithms. 
However, even with an optimum choice of depth 
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Figure 6: Influence of the regularization parameter 
(radial test-dipole at 71% relative eccentricity). 

normalization, MNLS underestimates the depth of 
sources at most eccentricities. The Laplacian 
smoothness approach (LORETA) at eccentricities 
above 70% over-estimates the source depths, which 
is due to the smoothness constraint, that cannot be 
well fulfilled at the source space boundaries. 
Nevertheless, the modified LORETA method shows 
the smallest overall localization errors and the best 
spatial resolution for the simulated data. 
Furthermore, due to the more indirect effect of the 
source strengths in the model term, lead-field 
normalization is less critical as compared to the 
standard MNLS algorithm. 

Compared to the linear methods, the non-linear 
approaches exhibit better localization accuracy and 
improved spatial resolution at high SNRs. The 
interpretation of the results is more difficult, since 
especially at small SNRs and with the L r sparse 
approach a few, relatively strong dipoles are 
reconstructed close to or around the true source 
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location and no smooth distributions with approxi¬ 
mately matching center positions can be found. 

For the more realistic low SNR-cases spatial 
resolution and localization accuracy decrease for all 
methods studied and at deep source positions the 
reconstructed current densities develop irregular 
shapes. The non-linear L r norm algorithm suffers 
most from low SNRs by showing scattered and 
spurious results for deep source locations. 

Adjusting the regularization parameter, that weights 
between data and model term, can be performed 
automatically for the MNLS and L r sparse 
approaches, if the noise of the data is known and 
specified correctly. All methods are found to be very 
sensitive to over-fitting the data in the realistic SNR 
cases, since then the algorithms will start to fit the 
noise or correlated background features of the data. 
Thus fixed or wrong choices of the regularization 
parameter will lead to spurious results or ghost 
sources, that explain more of the data than their 
signal contents. Small deviations between measured 
and calculated fields can always be achieved by 
adjusting small values of the regularization 
parameter. For current density methods, they are 
thus no longer a good criterion for the quality of the 
reconstruction result. It is essential to have a good 
estimate for the noise in the data and adjust the 
regularization parameter accordingly. 
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